Balken FEM mit Python#
In diesem Notebook wird die Lösung eines Balkenproblems mit der Methode der finiten Elemente (FEM) in Python implementiert. Die Lösung wird mit der analytischen Lösung verglichen und die Konvergenz der FEM-Lösung wird untersucht.
Die Differentialgleichung für die Durchbiegung \(w\) eines Balkens mit konstanter Biegesteifigkeit \(EI\) lautet:
Die schwache Form der Differentialgleichung wurde in der Vorlesung besprochen und lautet:
Hierbei stellen \(F_I\) und \(M_I\) die an den Knoten angreifenden Einzellasten und Momente dar
.
Das Balkenelement hat folgende interne Freiheitsgrade und angreifende Lasten:
Ziel ist es, die schwache Form der Differentialgleichung in ein algebraisches Gleichungssystem umzuwandeln:
Hierbei hat der Vektor \(\boldsymbol{F}\) der äußeren angreifenden Lasten nur Einträge an den Stellen, die mit den Durchbiegungen verknüpft sind. Der Vektor \(\boldsymbol{M}\) hat hingegen nur Einträge an den Stellen, die mit den Neigungen verknüpft sind.
Verschiebungsansatz#
Analog zum Stabelement führen wir die normierte Koordinate \(\xi=\frac{x}{\ell_e}\) ein. Die Ansatzfunktionen für die Durchbiegung \(w_h\) wird dann als Linearkombination der Formfunktionen \(N_I\) und der Freiheitsgrade \(\hat{w}_I\) geschrieben: $\( w_h = \sum_{I=1}^{6} N_I \hat{w}_I \)$
Als Formfunktionen \(N_I\) werden die Hermite-Polynome verwendet. Diese lauten:
import numpy as np
import plotly.graph_objects as go
# Define the x values
x = np.linspace(0, 1, 100)
# Define the shape functions - Hermite polynomials 4th order
N1 = 1 - 3*x**2 + 2*x**3
N2 = x - 2*x**2 + x**3
N3 = 3*x**2 - 2*x**3
N4 = -x**2 + x**3
# Create the figure
fig = go.Figure()
# Add traces for each shape function
fig.add_trace(go.Scatter(x=x, y=N1, mode='lines', name='N_1'))
fig.add_trace(go.Scatter(x=x, y=N2, mode='lines', name='N_2'))
fig.add_trace(go.Scatter(x=x, y=N3, mode='lines', name='N_3'))
fig.add_trace(go.Scatter(x=x, y=N4, mode='lines', name='N_4'))
fig.update_layout(
xaxis_title=r"$\xi$",
yaxis_title=r'$N(\xi)$',
legend=dict(x=1.05, y=1),
font=dict(family="Serif", size=15),
template='plotly_white',
xaxis=dict(showgrid=True, gridcolor='grey', gridwidth=0.6, griddash='dash'),
yaxis=dict(showgrid=True, gridcolor='grey', gridwidth=0.6, griddash='dash')
)
Element-Steifigkeitsmatrix für ein Balkenelemt#
Die Formfunktionen sind bzgl. der lokalen Koordinate \(\xi\) definiert. In der schwachen Formulierung wird die Ableitung der Formfunktionen bzgl. der globalen Koordinate \(x\) benötigt. Die Ableitung der Formfunktionen bzgl. der globalen Koordinate \(x\) ist gegeben durch: $\( \frac{\text{d} \xi}{\text{d} x} = \frac{1}{\ell_e} \qquad \rightarrow \qquad \text{d} x= \ell_e \text{d} \xi \)$
Wenn wir uns nun den ersten Term der schwachen Form herausgreifen und die Ableitungen bilden lauten diese:
Die Steifigkeitsmatrix \(\boldsymbol{K}_{IJ}\) wollen wir jetzt interaktiv berechnen.
Definition der Formfunktionen#
import sympy as sp
# Define the variables
xi,ell = sp.symbols('xi, ell')
# Define the Hermite polynomials
N1 = 1 - 3*xi**2 + 2*xi**3
N2 = (xi - 2*xi**2 + xi**3)*ell
N3 = 3*xi**2 - 2*xi**3
N4 = (-xi**2 + xi**3)*ell
N = sp.Matrix([[N1], [N2], [N3], [N4]])
N
Berechnung der Ableitungen#
Siehe schwache Form: $$
q_h(\xi) = (1-\xi) \hat{q}_1 + \xi \hat{q}_2 = \begin{pmatrix} 1-\xi & \xi \end{pmatrix} \begin{pmatrix} \hat{q}_1 \ \hat{q}_2 \end{pmatrix} $$
Dies wird jetzt in die schwache Form eingesetzt:
Den Vektor der Streckenlasten \(\boldsymbol{F}_Q\) berechnen wir nun interaktiv.
Ansatz für die Streckenlast#
q1,q2 = sp.symbols('q1 q2')
qh = sp.Matrix([q1,q2])
Nq1 = 1-xi
Nq2 = xi
Nq = sp.Matrix([Nq1, Nq2]).T
Nq*qh
Bildung des Integranten#
integrantq = N*Nq*qh*ell
integrantq
Integrieren der Funktion#
FQ = sp.integrate(integrantq,(xi,0,1))
sp.simplify(FQ)
Beispiel#
Für das dargestellte Beispiel eines statisch bestimmten Balkens mit konstanter Streckenlast lässt sich die Lösung der Verschiebung durch vierfache Integration der Differentialgleichung bestimmen. Die Lösung lautet:
Auch die Schnittgrößen können analytisch einfach bestimmt werden:
import sympy as sp
from sympy import Rational
from sympy import Eq
from sympy.plotting import plot, PlotGrid
x,EI,ell,q0,C_1,C_2,C_3,C_4 = sp.symbols('x,EI,ell,q0,C_1,C_2,C_3,C_4')
# Biegelinie
w = 1/EI*(Rational(1,24)*q0*x**4 +Rational(1,6)*C_1*x**3+Rational(1,2)*C_2*x**2+C_3*x+C_4)
# Bestimmung der Integrationskonstanten
eq1 = Eq(w.subs(x,0),0)
eq2 = Eq(w.subs(x,ell),0)
eq3 = Eq(-w.diff(x,2).subs(x,0),0)
eq4 = Eq(-w.diff(x,2).subs(x,ell),0)
solution = sp.solve((eq1,eq2,eq3,eq4),(C_1,C_2,C_3,C_4))
wfun = w.subs(solution)
Mfun = -wfun.diff(x,2)
Qfun = -wfun.diff(x,3)
# Plot der Biegelinie und Schnittgrößen
wp = w.subs(solution).subs({"ell":1,"q0":1,"EI":1})
M = -wp.diff(x,2)
Q = -wp.diff(x,3)
p1= plot(wp,(x,0,1),title=r'Durchbiegung $\frac{EI}{q_0} w$',show=False,xlabel=r'$\frac{x}{\ell}$',ylabel=r'$\frac{EI}{q_0} w$',line_color='black')
p2= plot(M,(x,0,1),title=r'Biegemoment $\frac{1}{q_0} M$',show=False,xlabel=r'$\frac{x}{\ell}$',ylabel=r'$\frac{1}{q_0}M$',line_color='red')
p3= plot(Q,(x,0,1),title=r'Querkraft $\frac{1}{q_0}Q$',show=False,xlabel=r'$\frac{x}{\ell}$',ylabel=r'$\frac{1}{q_0}Q$',line_color='blue')
PlotGrid(1,3,p1,p2,p3,size=(15,5));
Lösung mit der FEM#
import numpy as np
class BalkenFEM:
def __init__(self,numnp=2,numel=1):
self.numnp = numnp
self.numel = numel
self.dim = 2
self.nel = 2
self.numdof = self.numnp * self.dim
self.elements = np.zeros((self.numel,2),dtype=int)
self.eData = []
self.coords = np.zeros((self.numnp,2))
self.dof = np.zeros((self.numnp,self.dim))
self.eqind = np.zeros((self.numnp,self.dim))
self.Kges = np.zeros((self.numdof,self.numdof))
self.Fges = np.zeros((self.numnp,self.dim))
self.qL = np.zeros((self.numnp,1))
self.u = np.zeros((self.numnp,self.dim))
self.neumannBC = None
self.dirichletBC = None
self.Kuu = self.Kud = self.Kdu = self.Kdd = None
def setNodalCoords(self,coords):
self.coords[:,:] = coords[:,:]
def setElements(self,elements):
for i in range(self.numel):
for j in range(self.nel):
self.elements[i,j] = int(elements[i,j])
def setElementData(self,youngsmoduli,sma):
for i in range(self.numel):
self.eData.append({"sma":sma[i],"youngsmodulus":youngsmoduli[i]})
def setDirichletBoundaryCondition(self,dirichletNodes,dirichletDir,dirichletVal):
eqID = 1
for nodeID in range(self.numnp):
for dirID in range(self.dim):
self.eqind[nodeID,dirID] = eqID
eqID +=1
# Knoten an denen Dirichlet Randbedingungen vorgegeben werden
# müssen im Gleichungssystem nicht berücksichtigt werden.
# Deshalb erhalten sie negative Gleichungsnummern zur Identifikation
self.dirichletBC = np.zeros((len(dirichletNodes),3))
for i,(nodeID,dirId,dVal) in enumerate(zip(dirichletNodes,dirichletDir,dirichletVal)):
self.eqind[nodeID,dirId] = -self.eqind[nodeID,dirId]
self.dirichletBC[i,0] = nodeID
self.dirichletBC[i,1] = dirID
self.dirichletBC[i,2] = dVal
self.dof[nodeID,dirId] = dVal
def setExternalForces(self,neumannNodes,neumannDir,neumannVal):
for nodeID,nDir,nVal in zip(neumannNodes,neumannDir,neumannVal):
self.Fges[nodeID,nDir] = nVal
def setDistributedLoads(self,qloads):
self.qL[:] = qloads[:]
def _getElementDirector(self,elID):
nodeID1 = self.elements[elID,0]
nodeID2 = self.elements[elID,1]
director = self.coords[nodeID2,:] - self.coords[nodeID1,:]
return director
def _getElementstiffness_(self,elID):
director = self._getElementDirector(elID)
le = np.linalg.norm(director)
Ke = self.eData[elID]["youngsmodulus"]*self.eData[elID]["sma"]/(le**3) * np.ones((4,4))
Ke[0,0] *= 12
Ke[0,1] *= 6
Ke[0,2] *= -12
Ke[0,3] *= 6
Ke[1,0] *= 6
Ke[1,1] *= 4
Ke[1,2] *= -6
Ke[1,3] *= 2
Ke[2,0] *= -12
Ke[2,1] *= -6
Ke[2,2] *= 12
Ke[2,3] *= -6
Ke[3,0] *= 6
Ke[3,1] *= 2
Ke[3,2] *= -6
Ke[3,3] *= 4
Fq = np.zeros((4,1))
node1 = self.elements[elID,0]
node2 = self.elements[elID,1]
q1 = self.qL[node1]
q2 = self.qL[node2]
Fq[0] = le/20 * (7*q1+3*q2)
Fq[1] = le * (q1/20+q2/30)
Fq[2] = le/20 * (3*q1+7*q2)
Fq[3] = le * (-q1/30-q2/20)
return Ke,Fq
def secondDerivative(self,xi):
return np.array([6*(2*xi-1),2*(3*xi-2),6*(1-2*xi),2*(3*xi-1) ])
def thirdDerivative(self,xi):
return np.array([6*(2),2*(3),6*(-2),2*(3) ])
def computeMoment(self,n=10):
X = np.zeros(n*self.numel)
Mges = np.zeros(n*self.numel)
for elID in range(self.numel):
director = self._getElementDirector(elID)
le = np.linalg.norm(director)
node1 = self.elements[elID,0]
node2 =self.elements[elID,1]
x1 = self.coords[node1,0]
x2 = self.coords[node2,0]
X[elID*n:(elID+1)*n] = xlin = np.linspace(x1,x2,n)
xilin = np.linspace(0,1,n)
dofe = np.array([self.dof[node1,0],self.dof[node1,1],self.dof[node2,0],self.dof[node2,1]])
for i,xi in enumerate(xilin):
B=1.0/(le**2) * self.secondDerivative(xi)
# M = self.eData[elID]["youngsmodulus"]*self.eData[elID]["sma"] * B @ dofe.T
M = B @ dofe.T
Mges[elID*n+i] = -M
return X,Mges
def computeQuerkraft(self,n=10):
X = np.zeros(n*self.numel)
Qges = np.zeros(n*self.numel)
for elID in range(self.numel):
director = self._getElementDirector(elID)
le = np.linalg.norm(director)
node1 = self.elements[elID,0]
node2 =self.elements[elID,1]
x1 = self.coords[node1,0]
x2 = self.coords[node2,0]
X[elID*n:(elID+1)*n] = xlin = np.linspace(x1,x2,n)
xilin = np.linspace(0,1,n)
dofe = np.array([self.dof[node1,0],self.dof[node1,1],self.dof[node2,0],self.dof[node2,1]])
for i,xi in enumerate(xilin):
B=1.0/(le**3) * self.thirdDerivative(xi)
# M = self.eData[elID]["youngsmodulus"]*self.eData[elID]["sma"] * B @ dofe.T
Q = B @ dofe.T
Qges[elID*n+i] = -Q
return X,Qges
def getDeflection(self,x):
for elID in range(self.numel):
node1 = self.elements[elID,0]
node2 =self.elements[elID,1]
x1 = self.coords[node1,0]
x2 = self.coords[node2,0]
if (x-x1 > -1.e-12) and (x-x2 < 0.0 ):
director = self._getElementDirector(elID)
le = np.linalg.norm(director)
dofe = np.array([self.dof[node1,0],self.dof[node1,1],self.dof[node2,0],self.dof[node2,1]])
xi = (x-x1)/(x2-x1)
N = np.array([2*xi**3 - 3*xi**2 + 1, (xi**3 - 2*xi**2 + xi), -2*xi**3 + 3*xi**2, (xi**3 - xi**2)])
return N @ dofe.T
def assembleGlobalMatrix2D(self):
for elID in range(self.numel):
Ke,Fq = self._getElementstiffness_(elID)
for nodeI in range(self.nel):
globalNodeI = self.elements[elID,nodeI]
for dimI in range(self.dim):
rowID = int(np.sign(self.eqind[globalNodeI,dimI])*self.eqind[globalNodeI,dimI]-1)
self.Fges[globalNodeI,dimI] += Fq[(nodeI*self.dim)+dimI,0]
for nodeJ in range(self.nel):
globalNodeJ = self.elements[elID,nodeJ]
for dimJ in range(self.dim):
colID = int(np.sign(self.eqind[globalNodeJ,dimJ])*self.eqind[globalNodeJ,dimJ]-1)
self.Kges[rowID,colID] += Ke[(nodeI*self.dim)+dimI,(nodeJ*self.dim)+dimJ]
def solveSystem(self):
# Erzeuge Hiflsmatrizen um die Freiheitsgrade zu sortieren
numcdof = self.dirichletBC.shape[0]
numfdof = self.numdof - numcdof
self.Kuu = np.zeros((numfdof,numfdof))
self.Kud = np.zeros((numfdof,numcdof))
self.Kdd = np.zeros((numcdof,numcdof))
Fu = np.zeros((numfdof,1))
Fc = np.zeros((numcdof,1))
uc = np.zeros((numcdof,1))
eqf_inv = np.zeros((numfdof,2),dtype=int)
eqc_inv = np.zeros((numcdof,2),dtype=int)
Fvec = np.zeros((2*self.numnp,1))
# Sortiere Freiheitsgrade
icon = -1
ifree = -1
iIsFree = False
for nodeI in range(self.numnp):
for dirI in range(self.dim):
jcon = -1
jfree = -1
jIsFree = False
if self.eqind[nodeI,dirI]>=0:
iIsFree = True
ifree+=1
else:
iIsFree = False
icon +=1
eqI = int(np.sign(self.eqind[nodeI,dirI])*self.eqind[nodeI,dirI]-1)
Fvec[eqI] = self.Fges[nodeI,dirI]
if iIsFree:
eqf_inv[ifree,0] = nodeI
eqf_inv[ifree,1] = dirI
Fu[ifree]=self.Fges[nodeI,dirI]
else:
eqc_inv[icon,0] = nodeI
eqc_inv[icon,1] = dirI
Fc[icon] = self.Fges[nodeI,dirI]
uc[icon] = self.dof[nodeI,dirI]
for nodeJ in range(self.numnp):
for dirJ in range(self.dim):
if self.eqind[nodeJ,dirJ] >= 0 :
jIsFree = True
jfree +=1
else:
jIsFree = False
jcon += 1
eqJ = int(np.sign(self.eqind[nodeJ,dirJ])*self.eqind[nodeJ,dirJ]-1)
if iIsFree and jIsFree:
self.Kuu[ifree,jfree] = self.Kges[eqI,eqJ]
if iIsFree and not jIsFree:
self.Kud[ifree,jcon] = self.Kges[eqI,eqJ]
if not iIsFree and not jIsFree:
self.Kdd[icon,jcon] = self.Kges[eqI,eqJ]
# display("Fvec: ",Fvec)
# Löse Gleichungssystem
# K_uu*u = F-K_ud*u_d
rhs = Fu-self.Kud @ uc
# display("RHS: ",rhs)
# display("Kuu:",self.Kuu)
u = np.linalg.solve(self.Kuu,rhs)
ieq = 0
# Kopiere zu globalen Freiheitsgraden
for nodeI in range(numfdof):
self.dof[eqf_inv[nodeI,0],eqf_inv[nodeI,1]] = u[ieq,0]
ieq += 1
# Bestimme Reaktionskräfte
# R = K_du * u + K_dd*u_c - F_c
rF = self.Kud.T@u+self.Kdd@uc -Fc
# Kopiere zu den globalen Kräften
ieq = 0
for nodeI in range(numcdof):
self.Fges[eqc_inv[nodeI,0],eqc_inv[nodeI,1]] = rF[ieq,0]
ieq += 1
Set Up Problem#
numnp = 3 # Anzahl der Knoten im System
numel = numnp -1 # Anzahl der Elemente im System
ell = 8000 # Länge des Balkens in mm
# Erzeugen des FE-Objektes
Lastfall = BalkenFEM(numnp=numnp,numel=numel)
# Definiere die Koordinaten der beiden Randknoten
X1 = np.array([0.0,0.0]) #mm
X2 = np.array([ell,0.0]) #mm
# Generate interpolation points
t = np.linspace(0, 1, numnp)
X = X1 + t[:, np.newaxis] * (X2 - X1)
# print(f"Coordinates of the nodes: \n{X}")
Lastfall.setNodalCoords(X)
# Generate Elements
IX = np.column_stack((np.arange(0, numnp-1), np.arange(1, numnp)))
# print(f"Connectivity table: \n{IX}")
Lastfall.setElements(IX)
# Material properties and cross section
Emod = 210000 # MPa
sma = 77.67*10**4 # mm^4
Ilist = [sma for i in range(numel)]
Elist = [Emod for i in range(numel)]
Lastfall.setElementData(Elist, Ilist)
Definiere Randbedingungen#
# Dirichlet Randbedingungen (Verschiebung)
dnodes = [0, numnp-1]
ddir = [0,0]
dval = [0,0]
Lastfall.setDirichletBoundaryCondition(dnodes,ddir,dval)
# Neumann Randbedingungen (Kraft)
fnodes = [i for i in range(numnp)] # forces
mnodes = [i for i in range(numnp)] # moments
fdir1 = [0 for i in range(numnp)]
mdir2 = [1 for i in range(numnp)]
fval1 = [0 for i in range(numnp)]
mval2 = [0 for i in range(numnp)]
q0 = -10 # N/mm
# fval1[2] = q0*ell
Lastfall.setExternalForces(fnodes+mnodes,fdir1+mdir2,fval1+mval2)
# Verteilte Lasten
q=q0*np.ones((numnp,1))
Lastfall.setDistributedLoads(q)
Plot Problem#
import matplotlib.pyplot as plt
import scienceplots
plt.style.use(['science','grid'])
wanalytisch = sp.lambdify(x,wfun.subs({"EI":Emod*sma,"q0":q0,"ell":ell}),"numpy")
def plotBalken(X,q,fval1,u=None):
fix,ax = plt.subplots(figsize=(10,5))
if u is not None:
ax.plot(X[:,0],X[:,1]+u,c='b')
xline = np.linspace(0,ell,100)
ax.plot(xline,wanalytisch(xline),'b--')
ax.legend(["w - FEM","w - analytisch"])
ax.scatter(X[:,0],X[:,1],c='k',s=100)
ax.plot(X[:,0],X[:,1],c='k')
ax.quiver(X[:,0], X[:,1], [0 for i in range(numnp)],q,color='red')
ax.quiver(X[:,0], X[:,1], [0 for i in range(numnp)],fval1,color='red')
ax.set_xlabel('x [mm]')
ax.set_ylabel('w [mm]')
ax.set_title('Balken')
ax.grid(True)
plt.axis('equal')
plt.show()
plotBalken(X,q,fval1)
/home/steffen/anaconda3/envs/rise/lib/python3.10/site-packages/matplotlib/quiver.py:695: RuntimeWarning:
divide by zero encountered in scalar divide
/home/steffen/anaconda3/envs/rise/lib/python3.10/site-packages/matplotlib/quiver.py:695: RuntimeWarning:
invalid value encountered in multiply
Solve Problem#
Lastfall.assembleGlobalMatrix2D()
Lastfall.solveSystem()
plotBalken(X,q,fval1,Lastfall.dof[:,0])
wa = wanalytisch(4000)
wFEM = Lastfall.getDeflection(4000)
error = (wa-wFEM)/wa*100
print(f"Durchbiegung analytisch: {wa}")
print(f"Durchbiegung FEM: {wFEM}")
print(f"Fehler: {error}%")
Durchbiegung analytisch: -3269.8371825448057
Durchbiegung FEM: -3269.837182544792
Fehler: 4.172201784058398e-13%
Bestimmung der Schnittlasten#
Schnittmoment#
xm,M = Lastfall.computeMoment(n=10)
Manalytisch = sp.lambdify(x,Mfun.subs({"EI":Emod*sma,"q0":q0,"ell":ell}),"numpy")
def plotMoment(x,M):
fix,ax = plt.subplots(figsize=(10,5))
ax.plot(x,M,'r')
ax.plot(x,Manalytisch(x),'r--')
ax.legend(["M - FEM","M - analytisch"])
plotMoment(xm,M)
Ma = Manalytisch(4000)
MFEM = np.min(M)
error = (Ma-MFEM)/Ma*100
print(f"Biegemoment analytisch: {Ma}")
print(f"Biegemoment FEM: {MFEM}")
print(f"Fehler: {error}%")
Biegemoment analytisch: -0.0004904755773817184
Biegemoment FEM: -0.0005722215069453389
Fehler: -16.666666666666835%
Schnittgröße: Querkraft#
xm,Q = Lastfall.computeQuerkraft(n=10)
Qanalytisch = sp.lambdify(x,Qfun.subs({"EI":Emod*sma,"q0":q0,"ell":ell}),"numpy")
def plotQuerkraft(x,Q):
fix,ax = plt.subplots(figsize=(10,5))
ax.plot(x,Q,'b')
ax.plot(x,Qanalytisch(x),'b--')
ax.legend(["Q - FEM","Q - analytisch"])
plotQuerkraft(xm,Q)
Qa = Qanalytisch(0)
QFEM = np.min(Q)
error = (Qa-QFEM)/Qa*100
print(f"Querkraft analytisch: {Qa}")
print(f"Querkraft FEM: {QFEM}")
print(f"Fehler: {error}%")
Querkraft analytisch: -2.45237788690859e-07
Querkraft FEM: -1.2261889434542992e-07
Fehler: 49.99999999999983%
Konvergenzstudie#
def setUpProblem(numnp):
numel = numnp -1
ell = 8000
Lastfall = BalkenFEM(numnp=numnp,numel=numel)
X1 = np.array([0.0,0.0]) #mm
X2 = np.array([ell,0.0]) #mm
# Generate interpolation points
t = np.linspace(0, 1, numnp)
X = X1 + t[:, np.newaxis] * (X2 - X1)
# print(f"Coordinates of the nodes: \n{X}")
Lastfall.setNodalCoords(X)
# Generate Elements
IX = np.column_stack((np.arange(0, numnp-1), np.arange(1, numnp)))
# print(f"Connectivity table: \n{IX}")
Lastfall.setElements(IX)
# Material properties and cross section
Emod = 210000 # MPa
sma = 77.67*10**4 # mm^4
Ilist = [sma for i in range(numel)]
Elist = [Emod for i in range(numel)]
Lastfall.setElementData(Elist, Ilist)
# Dirichlet Randbedingungen (Verschiebung)
dnodes = [0, numnp-1]
ddir = [0,0]
dval = [0,0]
Lastfall.setDirichletBoundaryCondition(dnodes,ddir,dval)
# Neumann Randbedingungen (Kraft)
fnodes = [i for i in range(numnp)] # forces
mnodes = [i for i in range(numnp)] # moments
fdir1 = [0 for i in range(numnp)]
mdir2 = [1 for i in range(numnp)]
fval1 = [0 for i in range(numnp)]
mval2 = [0 for i in range(numnp)]
q0 = -10 # N/mm
# fval1[2] = q0*ell
Lastfall.setExternalForces(fnodes+mnodes,fdir1+mdir2,fval1+mval2)
# Verteilte Lasten
q=q0*np.ones((numnp,1))
Lastfall.setDistributedLoads(q)
Lastfall.assembleGlobalMatrix2D()
Lastfall.solveSystem()
xm,Q = Lastfall.computeQuerkraft(n=3)
xm,M = Lastfall.computeMoment(n=3)
return np.min(Q),np.min(M),Lastfall
npoints = np.linspace(2,202,51)
QFEM = np.zeros((51,))
MFEM = np.zeros((51,))
ErrorM= np.zeros((51,))
ErrorQ= np.zeros((51,))
for i,p in enumerate(npoints):
QQ,MM,LC=setUpProblem(int(p))
QFEM[i]=QQ
MFEM[i]=MM
ErrorQ[i] = (Qa-QFEM[i])/Qa*100
ErrorM[i] = (Ma-MFEM[i])/Ma*100
import pandas as pd
df = pd.DataFrame({"M":MFEM,"error M":ErrorM,"Q":QFEM,"error Q":ErrorQ},index=npoints)
df
| M | error M | Q | error Q | |
|---|---|---|---|---|
| 2.0 | -0.000327 | 33.333333 | -0.000000e+00 | 100.000000 |
| 6.0 | -0.000484 | 1.333333 | -1.961902e-07 | 20.000000 |
| 10.0 | -0.000488 | 0.411523 | -2.179891e-07 | 11.111111 |
| 14.0 | -0.000490 | 0.197239 | -2.263733e-07 | 7.692308 |
| 18.0 | -0.000490 | 0.115340 | -2.308120e-07 | 5.882353 |
| 22.0 | -0.000490 | 0.075586 | -2.335598e-07 | 4.761905 |
| 26.0 | -0.000490 | 0.053333 | -2.354283e-07 | 4.000000 |
| 30.0 | -0.000490 | 0.039635 | -2.367813e-07 | 3.448276 |
| 34.0 | -0.000490 | 0.030609 | -2.378063e-07 | 3.030303 |
| 38.0 | -0.000490 | 0.024349 | -2.386097e-07 | 2.702703 |
| 42.0 | -0.000490 | 0.019829 | -2.392564e-07 | 2.439024 |
| 46.0 | -0.000490 | 0.016461 | -2.397881e-07 | 2.222222 |
| 50.0 | -0.000490 | 0.013883 | -2.402329e-07 | 2.040816 |
| 54.0 | -0.000490 | 0.011867 | -2.406107e-07 | 1.886792 |
| 58.0 | -0.000490 | 0.010260 | -2.409354e-07 | 1.754386 |
| 62.0 | -0.000490 | 0.008958 | -2.412175e-07 | 1.639344 |
| 66.0 | -0.000490 | 0.007890 | -2.414649e-07 | 1.538462 |
| 70.0 | -0.000490 | 0.007001 | -2.416836e-07 | 1.449275 |
| 74.0 | -0.000490 | 0.006255 | -2.418784e-07 | 1.369863 |
| 78.0 | -0.000490 | 0.005622 | -2.420529e-07 | 1.298701 |
| 82.0 | -0.000490 | 0.005081 | -2.422102e-07 | 1.234568 |
| 86.0 | -0.000490 | 0.004614 | -2.423526e-07 | 1.176471 |
| 90.0 | -0.000490 | 0.004208 | -2.424823e-07 | 1.123595 |
| 94.0 | -0.000490 | 0.003854 | -2.426008e-07 | 1.075269 |
| 98.0 | -0.000490 | 0.003543 | -2.427096e-07 | 1.030928 |
| 102.0 | -0.000490 | 0.003268 | -2.428097e-07 | 0.990099 |
| 106.0 | -0.000490 | 0.003023 | -2.429022e-07 | 0.952381 |
| 110.0 | -0.000490 | 0.002806 | -2.429879e-07 | 0.917431 |
| 114.0 | -0.000490 | 0.002611 | -2.430675e-07 | 0.884956 |
| 118.0 | -0.000490 | 0.002435 | -2.431417e-07 | 0.854701 |
| 122.0 | -0.000490 | 0.002277 | -2.432110e-07 | 0.826446 |
| 126.0 | -0.000490 | 0.002133 | -2.432759e-07 | 0.800000 |
| 130.0 | -0.000490 | 0.002003 | -2.433367e-07 | 0.775194 |
| 134.0 | -0.000490 | 0.001884 | -2.433939e-07 | 0.751880 |
| 138.0 | -0.000490 | 0.001776 | -2.434477e-07 | 0.729927 |
| 142.0 | -0.000490 | 0.001677 | -2.434985e-07 | 0.709220 |
| 146.0 | -0.000490 | 0.001586 | -2.435465e-07 | 0.689655 |
| 150.0 | -0.000490 | 0.001501 | -2.435919e-07 | 0.671141 |
| 154.0 | -0.000490 | 0.001424 | -2.436349e-07 | 0.653595 |
| 158.0 | -0.000490 | 0.001352 | -2.436758e-07 | 0.636943 |
| 162.0 | -0.000490 | 0.001286 | -2.437146e-07 | 0.621118 |
| 166.0 | -0.000490 | 0.001225 | -2.437515e-07 | 0.606061 |
| 170.0 | -0.000490 | 0.001166 | -2.437867e-07 | 0.591715 |
| 174.0 | -0.000490 | 0.001113 | -2.438202e-07 | 0.578034 |
| 178.0 | -0.000490 | 0.001064 | -2.438523e-07 | 0.564972 |
| 182.0 | -0.000490 | 0.001016 | -2.438829e-07 | 0.552486 |
| 186.0 | -0.000490 | 0.000974 | -2.439122e-07 | 0.540540 |
| 190.0 | -0.000490 | 0.000935 | -2.439402e-07 | 0.529102 |
| 194.0 | -0.000490 | 0.000895 | -2.439671e-07 | 0.518135 |
| 198.0 | -0.000490 | 0.000859 | -2.439929e-07 | 0.507614 |
| 202.0 | -0.000490 | 0.000824 | -2.440177e-07 | 0.497512 |
fig, ax = plt.subplots(1,2,figsize=(12,5))
ax[0].plot(df.index,df["M"])
ax[0].set_title("Konvergenz Schnittmoment")
ax[0].set_xlabel("Anzahl Knoten")
ax[0].set_ylabel("M")
ax[1].plot(df.index,df["Q"])
ax[1].set_title("Konvergenz SchnittKraft Q")
ax[1].set_xlabel("Anzahl Knoten")
ax[1].set_ylabel("Q")
Text(0, 0.5, 'Q')
fig, ax = plt.subplots(1,2,figsize=(12,5))
ax[0].plot(df.index,df["error M"])
ax[0].set_title("Konvergenz Fehler des Schnittmoments")
ax[0].set_xlabel("Anzahl Knoten")
ax[0].set_ylabel("error M")
ax[0].set_yscale("log") # Set y-axis to logarithmic scale
ax[0].set_xscale("log") # Set x-axis to logarithmic scale
ax[1].plot(df.index,df["error Q"])
ax[1].set_title("Konvergenz Fehler SchnittKraft Q")
ax[1].set_xlabel("Anzahl Knoten")
ax[1].set_ylabel("error Q")
ax[1].set_yscale("log") # Set y-axis to logarithmic scale
ax[1].set_xscale("log") # Set x-axis to logarithmic scale
Analyse der Steifigkeitsmatrix#
matrix = LC.Kges
plt.figure(figsize=(6, 6))
plt.spy(matrix, markersize=10) # You can adjust the markersize as needed
plt.title('Sparsity Pattern of the Matrix')
plt.xlabel('Column Index')
plt.ylabel('Row Index')
plt.grid()
plt.show()